home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / C / Applications / Python 1.3.3 / Python 133 SRC / Modules / mpzmodule.c < prev    next >
Text File  |  1995-12-21  |  41KB  |  1,810 lines
  1. /***********************************************************
  2. Copyright 1991-1995 by Stichting Mathematisch Centrum, Amsterdam,
  3. The Netherlands.
  4.  
  5.                         All Rights Reserved
  6.  
  7. Permission to use, copy, modify, and distribute this software and its 
  8. documentation for any purpose and without fee is hereby granted, 
  9. provided that the above copyright notice appear in all copies and that
  10. both that copyright notice and this permission notice appear in 
  11. supporting documentation, and that the names of Stichting Mathematisch
  12. Centrum or CWI not be used in advertising or publicity pertaining to
  13. distribution of the software without specific, written prior permission.
  14.  
  15. STICHTING MATHEMATISCH CENTRUM DISCLAIMS ALL WARRANTIES WITH REGARD TO
  16. THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
  17. FITNESS, IN NO EVENT SHALL STICHTING MATHEMATISCH CENTRUM BE LIABLE
  18. FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  19. WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  20. ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT
  21. OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  22.  
  23. ******************************************************************/
  24.  
  25. /* MPZ module */
  26.  
  27. /* This module provides an interface to an alternate Multi-Precision
  28.    library, GNU MP in this case */
  29.  
  30. /* XXX note: everywhere where mpz_size is called,
  31.    sizeof (limb) == sizeof (long)  has been assumed. */
  32.    
  33.  
  34. /* MPZ objects */
  35.  
  36. #include "allobjects.h"
  37.  
  38. #include <assert.h>
  39. #include <sys/types.h>        /* For size_t */
  40.  
  41. /*
  42. **    These are the cpp-flags used in this file...
  43. **
  44. **
  45. ** MPZ_MDIV_BUG        works around the mpz_m{div,mod,...} routines.
  46. **             This bug has been fixed in a later release of
  47. **             GMP.
  48. ** 
  49. ** MPZ_GET_STR_BUG    mpz_get_str corrupts memory, seems to be fixed
  50. **             in a later release
  51. ** 
  52. ** MPZ_DEBUG        generates a bunch of diagnostic messages
  53. ** 
  54. ** MPZ_SPARE_MALLOC    if set, results in extra code that tries to
  55. **             minimize the creation of extra objects.
  56. ** 
  57. ** MPZ_TEST_DIV        extra diagnostic output on stderr, when division
  58. **             routines are involved
  59. ** 
  60. ** MPZ_LIB_DOES_CHECKING    if set, assumes that mpz library doesn't call
  61. **             alloca with arg < 0 (when casted to a signed
  62. **             integral type).
  63. ** 
  64. ** MPZ_CONVERSIONS_AS_METHODS    if set, presents the conversions as
  65. **             methods. e.g., `mpz(5).long() == 5L'
  66. **             Later, Guido provided an interface to the
  67. **             standard functions. So this flag has no been
  68. **             cleared, and `long(mpz(5)) == 5L'
  69. ** 
  70. ** MP_TEST_ALLOC    If set, you would discover why MPZ_GET_STR_BUG
  71. **            is needed
  72. ** 
  73. ** MAKEDUMMYINT        Must be set if dynamic linking will be used
  74. */
  75.  
  76.  
  77. /*
  78. ** IMHO, mpz_m{div,mod,divmod}() do the wrong things when the denominator < 0
  79. ** I assume that this will be fixed in a future release
  80. */
  81. /*#define MPZ_MDIV_BUG fixed the (for me) nexessary parts in libgmp.a */
  82. /*
  83. ** IMO, mpz_get_str() assumes a bit too large target space, if he doesn't
  84. ** allocate it himself
  85. */
  86. #define MPZ_GET_STR_BUG
  87.  
  88. #include "gmp.h"
  89. typedef struct {
  90.     OB_HEAD
  91.         MP_INT    mpz;        /* the actual number */
  92. } mpzobject;
  93.  
  94. staticforward typeobject MPZtype;
  95.  
  96. #define is_mpzobject(v)        ((v)->ob_type == &MPZtype)
  97.  
  98. static const char initialiser_name[] = "mpz";
  99.  
  100. /* #define MPZ_DEBUG */
  101.  
  102. static mpzobject *
  103. newmpzobject()
  104. {
  105.     mpzobject *mpzp;
  106.  
  107.  
  108. #ifdef MPZ_DEBUG
  109.     fputs( "mpz_object() called...\n", stderr );
  110. #endif /* def MPZ_DEBUG */
  111.     mpzp = NEWOBJ(mpzobject, &MPZtype);
  112.     if (mpzp == NULL)
  113.         return NULL;
  114.  
  115.     mpz_init(&mpzp->mpz);    /* actual initialisation */
  116.     return mpzp;
  117. } /* newmpzobject() */
  118.  
  119. #ifdef MPZ_GET_STR_BUG
  120. #include "gmp-impl.h"
  121. #include "longlong.h"
  122. #endif /* def MPZ_GET_STR_BUG */
  123.  
  124. static object *
  125. mpz_format(objp, base, withname)
  126.     object *objp;
  127.     int base;
  128.     unsigned char withname;
  129. {
  130.     mpzobject *mpzp = (mpzobject *)objp;
  131.     stringobject *strobjp;
  132.     int i;
  133.     int cmpres;
  134.     int taglong;
  135.     char *cp;
  136.     char prefix[5], *tcp;
  137.  
  138.  
  139.     tcp = &prefix[0];
  140.  
  141.     if (mpzp == NULL || !is_mpzobject(mpzp)) {
  142.         err_badcall();
  143.         return NULL;
  144.     }
  145.  
  146.     assert(base >= 2 && base <= 36);
  147.  
  148.     if (withname)
  149.         i = strlen(initialiser_name) + 2; /* e.g. 'mpz(' + ')' */
  150.     else
  151.         i = 0;
  152.  
  153.     if ((cmpres = mpz_cmp_si(&mpzp->mpz, 0L)) == 0)
  154.         base = 10;    /* '0' in every base, right */
  155.     else if (cmpres < 0) {
  156.         *tcp++ = '-';
  157.         i += 1;        /* space to hold '-' */
  158.     }
  159.  
  160. #ifdef MPZ_DEBUG
  161.     fprintf(stderr, "mpz_format: mpz_sizeinbase %d\n",
  162.         (int)mpz_sizeinbase(&mpzp->mpz, base));
  163. #endif /* def MPZ_DEBUG */
  164. #ifdef MPZ_GET_STR_BUG
  165.     i += ((size_t) abs(mpzp->mpz.size) * BITS_PER_MP_LIMB
  166.           * __mp_bases[base].chars_per_bit_exactly) + 1;
  167. #else /* def MPZ_GET_STR_BUG */
  168.     i += (int)mpz_sizeinbase(&mpzp->mpz, base);
  169. #endif /* def MPZ_GET_STR_BUG else */
  170.  
  171.     if (base == 16) {
  172.         *tcp++ = '0';
  173.         *tcp++ = 'x';
  174.         i += 2;        /* space to hold '0x' */
  175.     }
  176.     else if (base == 8) {
  177.         *tcp++ = '0';
  178.         i += 1;        /* space to hold the extra '0' */
  179.     }
  180.     else if (base > 10) {
  181.         *tcp++ = '0' + base / 10;
  182.         *tcp++ = '0' + base % 10;
  183.         *tcp++ = '#';
  184.         i += 3;        /* space to hold e.g. '12#' */
  185.     }
  186.     else if (base < 10) {
  187.         *tcp++ = '0' + base;
  188.         *tcp++ = '#';
  189.         i += 2;        /* space to hold e.g. '6#' */
  190.     }
  191.  
  192.     /*
  193.     ** the following code looks if we need a 'L' attached to the number
  194.     ** it will also attach an 'L' to the value -0x80000000
  195.     */
  196.     taglong = 0;
  197.     if (mpz_size(&mpzp->mpz) > 1
  198.         || (long)mpz_get_ui(&mpzp->mpz) < 0L) {
  199.         taglong = 1;
  200.         i += 1;        /* space to hold 'L' */
  201.     }
  202.  
  203. #ifdef MPZ_DEBUG
  204.     fprintf(stderr, "mpz_format: requesting string size %d\n", i);
  205. #endif /* def MPZ_DEBUG */    
  206.     if ((strobjp = (stringobject *)newsizedstringobject((char *)0, i))
  207.         == NULL)
  208.         return NULL;
  209.  
  210.     /* get the beginning of the string memory and start copying things */
  211.     cp = GETSTRINGVALUE(strobjp);
  212.     if (withname) {
  213.         strcpy(cp, initialiser_name);
  214.         cp += strlen(initialiser_name);
  215.         *cp++ = '('; /*')'*/
  216.     }
  217.  
  218.     /* copy the already prepared prefix; e.g. sign and base indicator */
  219.     *tcp = '\0';
  220.     strcpy(cp, prefix);
  221.     cp += tcp - prefix;
  222.  
  223.     /* since' we have the sign already, let the lib think it's a positive
  224.         number */
  225.     if (cmpres < 0)
  226.         mpz_neg(&mpzp->mpz,&mpzp->mpz);    /* hack Hack HAck HACk HACK */
  227.     (void)mpz_get_str(cp, base, &mpzp->mpz);
  228.     if (cmpres < 0)
  229.         mpz_neg(&mpzp->mpz,&mpzp->mpz);    /* hack Hack HAck HACk HACK */
  230. #ifdef MPZ_DEBUG
  231.     fprintf(stderr, "mpz_format: base (ultim) %d, mpz_get_str: %s\n",
  232.         base, cp);
  233. #endif /* def MPZ_DEBUG */
  234.     cp += strlen(cp);
  235.  
  236.     if (taglong)
  237.         *cp++ = 'L';
  238.     if (withname)
  239.         *cp++ = /*'('*/ ')';
  240.  
  241.     *cp = '\0';
  242.  
  243. #ifdef MPZ_DEBUG
  244.     fprintf(stderr,
  245.         "mpz_format: cp (str end) 0x%x, begin 0x%x, diff %d, i %d\n",
  246.         cp, GETSTRINGVALUE(strobjp), cp - GETSTRINGVALUE(strobjp), i);
  247. #endif /* def MPZ_DEBUG */    
  248.     assert(cp - GETSTRINGVALUE(strobjp) <= i);
  249.  
  250.     if (cp - GETSTRINGVALUE(strobjp) != i) {
  251.         strobjp->ob_size -= i - (cp - GETSTRINGVALUE(strobjp));
  252.     }
  253.  
  254.     return (object *)strobjp;
  255. } /* mpz_format() */
  256.  
  257. /* MPZ methods */
  258.  
  259. static void
  260. mpz_dealloc(mpzp)
  261.     mpzobject *mpzp;
  262. {
  263. #ifdef MPZ_DEBUG
  264.     fputs( "mpz_dealloc() called...\n", stderr );
  265. #endif /* def MPZ_DEBUG */
  266.     mpz_clear(&mpzp->mpz);
  267.     DEL(mpzp);
  268. } /* mpz_dealloc() */
  269.  
  270.  
  271. /* pointers to frequently used values 0, 1 and -1 */
  272. static mpzobject *mpz_value_zero, *mpz_value_one, *mpz_value_mone;
  273.  
  274. static int
  275. mpz_compare(a, b)
  276.     mpzobject *a, *b;
  277. {
  278.     int cmpres;
  279.  
  280.  
  281.     /* guido sez it's better to return -1, 0 or 1 */
  282.     return (cmpres = mpz_cmp( &a->mpz, &b->mpz )) == 0 ? 0
  283.         : cmpres > 0 ? 1 : -1;
  284. } /* mpz_compare() */
  285.  
  286. static object *
  287. mpz_addition(a, b)
  288.     mpzobject *a;
  289.     mpzobject *b;
  290. {
  291.     mpzobject *z;
  292.  
  293.     
  294. #ifdef MPZ_SPARE_MALLOC
  295.     if (mpz_cmp_ui(&a->mpz, (unsigned long int)0) == 0) {
  296.         INCREF(b);
  297.         return (object *)b;
  298.     }
  299.  
  300.     if (mpz_cmp_ui(&b->mpz, (unsigned long int)0) == 0) {
  301.         INCREF(a);
  302.         return (object *)a;
  303.     }
  304. #endif /* def MPZ_SPARE_MALLOC */
  305.  
  306.     if ((z = newmpzobject()) == NULL)
  307.         return NULL;
  308.     
  309.     mpz_add(&z->mpz, &a->mpz, &b->mpz);
  310.     return (object *)z;
  311. } /* mpz_addition() */
  312.  
  313. static object *
  314. mpz_substract(a, b)
  315.     mpzobject *a;
  316.     mpzobject *b;
  317. {
  318.     mpzobject *z;
  319.  
  320.     
  321. #ifdef MPZ_SPARE_MALLOC
  322.     if (mpz_cmp_ui(&b->mpz, (unsigned long int)0) == 0) {
  323.         INCREF(a);
  324.         return (object *)a;
  325.     }
  326. #endif /* MPZ_SPARE_MALLOC */    
  327.  
  328.     if ((z = newmpzobject()) == NULL)
  329.         return NULL;
  330.  
  331.     mpz_sub(&z->mpz, &a->mpz, &b->mpz);
  332.     return (object *)z;
  333. } /* mpz_substract() */
  334.  
  335. static object *
  336. mpz_multiply(a, b)
  337.     mpzobject *a;
  338.     mpzobject *b;
  339. {
  340. #ifdef MPZ_SPARE_MALLOC
  341.     int cmpres;
  342. #endif /* def MPZ_SPARE_MALLOC */
  343.     mpzobject *z;
  344.  
  345.  
  346. #ifdef MPZ_SPARE_MALLOC
  347.     if ((cmpres = mpz_cmp_ui(&a->mpz, (unsigned long int)0)) == 0) {
  348.         INCREF(mpz_value_zero);
  349.         return (object *)mpz_value_zero;
  350.     }
  351.     if (cmpres > 0 && mpz_cmp_ui(&a->mpz, (unsigned long int)1) == 0) {
  352.         INCREF(b);
  353.         return (object *)b;
  354.     }
  355.  
  356.     if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long_int)0)) == 0) {
  357.         INCREF(mpz_value_zero);
  358.         return (object *)mpz_value_zero;
  359.     }
  360.     if (cmpres > 0 && mpz_cmp_ui(&b->mpz, (unsigned long int)1) == 0) {
  361.         INCREF(a);
  362.         return (object *)a;
  363.     }
  364. #endif /* MPZ_SPARE_MALLOC */
  365.  
  366.     if ((z = newmpzobject()) == NULL)
  367.         return NULL;
  368.  
  369.     mpz_mul( &z->mpz, &a->mpz, &b->mpz );
  370.     return (object *)z;
  371.     
  372. } /* mpz_multiply() */
  373.  
  374. static object *
  375. mpz_divide(a, b)
  376.     mpzobject *a;
  377.     mpzobject *b;
  378. {
  379. #ifdef MPZ_SPARE_MALLOC
  380.     int cmpres;
  381. #endif /* def MPZ_SPARE_MALLOC */
  382.     mpzobject *z;
  383.  
  384.  
  385.     if ((
  386. #ifdef MPZ_SPARE_MALLOC
  387.          cmpres =
  388. #endif /* def MPZ_SPARE_MALLOC */
  389.          mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
  390.         err_setstr(ZeroDivisionError, "mpz./ by zero");
  391.         return NULL;
  392.     }
  393. #ifdef MPZ_SPARE_MALLOC
  394.     if (cmpres > 0 && mpz_cmp_ui(&b->mpz(unsigned long int)1) == 0) {
  395.         INCREF(a);
  396.         return (object *)a;
  397.     }
  398. #endif /* def MPZ_SPARE_MALLOC */
  399.  
  400.     if ((z = newmpzobject()) == NULL)
  401.         return NULL;
  402.  
  403. #ifdef MPZ_TEST_DIV
  404.     fputs("mpz_divide:  div result", stderr);
  405.     mpz_div(&z->mpz, &a->mpz, &b->mpz);
  406.     mpz_out_str(stderr, 10, &z->mpz);
  407.     putc('\n', stderr);
  408. #endif /* def MPZ_TEST_DIV */
  409. #ifdef MPZ_MDIV_BUG
  410.     if ((mpz_cmp_ui(&a->mpz, (unsigned long int)0) < 0)
  411.         != (mpz_cmp_ui(&b->mpz, (unsigned long int)0) < 0)) {
  412.         /*
  413.         ** numerator has other sign than denominator: we have
  414.         ** to look at the remainder for a correction, since mpz_mdiv
  415.         ** also calls mpz_divmod, I can as well do it myself
  416.         */
  417.         MP_INT tmpmpz;
  418.  
  419.  
  420.         mpz_init(&tmpmpz);
  421.         mpz_divmod(&z->mpz, &tmpmpz, &a->mpz, &b->mpz);
  422.  
  423.         if (mpz_cmp_ui(&tmpmpz, (unsigned long int)0) != 0)
  424.             mpz_sub_ui(&z->mpz, &z->mpz, (unsigned long int)1);
  425.  
  426.         mpz_clear(&tmpmpz);
  427.     }
  428.     else
  429.         mpz_div(&z->mpz, &a->mpz, &b->mpz);
  430.         /* the ``naive'' implementation does it right for operands
  431.            having the same sign */
  432.  
  433. #else /* def MPZ_MDIV_BUG */
  434.     mpz_mdiv(&z->mpz, &a->mpz, &b->mpz);
  435. #endif /* def MPZ_MDIV_BUG else */
  436. #ifdef MPZ_TEST_DIV
  437.     fputs("mpz_divide: mdiv result", stderr);
  438.     mpz_out_str(stderr, 10, &z->mpz);
  439.     putc('\n', stderr);
  440. #endif /* def MPZ_TEST_DIV */
  441.     return (object *)z;
  442.     
  443. } /* mpz_divide() */
  444.  
  445. static object *
  446. mpz_remainder(a, b)
  447.     mpzobject *a;
  448.     mpzobject *b;
  449. {
  450. #ifdef MPZ_SPARE_MALLOC
  451.     int cmpres;
  452. #endif /* def MPZ_SPARE_MALLOC */    
  453.     mpzobject *z;
  454.  
  455.     
  456.     if ((
  457. #ifdef MPZ_SPARE_MALLOC         
  458.          cmpres =
  459. #endif /* def MPZ_SPARE_MALLOC */    
  460.          mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
  461.         err_setstr(ZeroDivisionError, "mpz.% by zero");
  462.         return NULL;
  463.     }
  464. #ifdef MPZ_SPARE_MALLOC
  465.     if (cmpres > 0) {
  466.         if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)2)) == 0) {
  467.             INCREF(mpz_value_one);
  468.             return (object *)mpz_value_one;
  469.         }
  470.         if (cmpres < 0) {
  471.             /* b must be 1 now */
  472.             INCREF(mpz_value_zero);
  473.             return (object *)mpz_value_zero;
  474.         }
  475.     }
  476. #endif /* def MPZ_SPARE_MALLOC */    
  477.  
  478.     if ((z = newmpzobject()) == NULL)
  479.         return NULL;
  480.  
  481. #ifdef MPZ_TEST_DIV
  482.     fputs("mpz_remain:  mod result", stderr);
  483.     mpz_mod(&z->mpz, &a->mpz, &b->mpz);
  484.     mpz_out_str(stderr, 10, &z->mpz);
  485.     putc('\n', stderr);
  486. #endif /* def MPZ_TEST_DIV */
  487. #ifdef MPZ_MDIV_BUG
  488.  
  489.     /* the ``naive'' implementation does it right for operands
  490.        having the same sign */
  491.     mpz_mod(&z->mpz, &a->mpz, &b->mpz);
  492.  
  493.     /* assumption: z, a and b all point to different locations */
  494.     if ((mpz_cmp_ui(&a->mpz, (unsigned long int)0) < 0)
  495.         != (mpz_cmp_ui(&b->mpz, (unsigned long int)0) < 0)
  496.         && mpz_cmp_ui(&z->mpz, (unsigned long int)0) != 0)
  497.         mpz_add(&z->mpz, &z->mpz, &b->mpz);
  498.         /*
  499.         ** numerator has other sign than denominator: we have
  500.         ** to look at the remainder for a correction, since mpz_mdiv
  501.         ** also calls mpz_divmod, I can as well do it myself
  502.         */
  503. #else /* def MPZ_MDIV_BUG */
  504.     mpz_mmod(&z->mpz, &a->mpz, &b->mpz);
  505. #endif /* def MPZ_MDIV_BUG else */
  506. #ifdef MPZ_TEST_DIV
  507.     fputs("mpz_remain: mmod result", stderr);
  508.     mpz_out_str(stderr, 10, &z->mpz);
  509.     putc('\n', stderr);
  510. #endif /* def MPZ_TEST_DIV */
  511.     return (object *)z;
  512.     
  513. } /* mpz_remainder() */
  514.  
  515. static object *
  516. mpz_div_and_mod(a, b)
  517.     mpzobject *a;
  518.     mpzobject *b;
  519. {
  520.     object *z = NULL;
  521.     mpzobject *x = NULL, *y = NULL;
  522.  
  523.  
  524.     if (mpz_cmp_ui(&b->mpz, (unsigned long int)0) == 0) {
  525.         err_setstr(ZeroDivisionError, "mpz.divmod by zero");
  526.         return NULL;
  527.     }
  528.  
  529.     if ((z = newtupleobject(2)) == NULL
  530.         || (x = newmpzobject()) == NULL
  531.         || (y = newmpzobject()) == NULL) {
  532.         XDECREF(z);
  533.         XDECREF(x);
  534.         XDECREF(y);
  535.         return NULL;
  536.     }
  537.  
  538. #ifdef MPZ_TEST_DIV
  539.     fputs("mpz_divmod:  dm  result", stderr);
  540.     mpz_divmod(&x->mpz, &y->mpz, &a->mpz, &b->mpz);
  541.     mpz_out_str(stderr, 10, &x->mpz);
  542.     putc('\n', stderr);
  543.     mpz_out_str(stderr, 10, &y->mpz);
  544.     putc('\n', stderr);
  545. #endif /* def MPZ_TEST_DIV */
  546. #ifdef MPZ_MDIV_BUG
  547.     mpz_divmod(&x->mpz, &y->mpz, &a->mpz, &b->mpz);
  548.     if ((mpz_cmp_ui(&a->mpz, (unsigned long int)0) < 0)
  549.         != (mpz_cmp_ui(&b->mpz, (unsigned long int)0) < 0)
  550.         && mpz_cmp_ui(&y->mpz, (unsigned long int)0) != 0) {
  551.         /*
  552.         ** numerator has other sign than denominator: we have
  553.         ** to look at the remainder for a correction.
  554.         */
  555.         mpz_add(&y->mpz, &y->mpz, &b->mpz);
  556.         mpz_sub_ui(&x->mpz, &x->mpz, (unsigned long int)1);
  557.     }
  558. #else /* def MPZ_MDIV_BUG */
  559.     mpz_mdivmod( &x->mpz, &y->mpz, &a->mpz, &b->mpz );
  560. #endif /* def MPZ_MDIV_BUG else */
  561. #ifdef MPZ_TEST_DIV
  562.     fputs("mpz_divmod: mdm  result", stderr);
  563.     mpz_out_str(stderr, 10, &x->mpz);
  564.     putc('\n', stderr);
  565.     mpz_out_str(stderr, 10, &y->mpz);
  566.     putc('\n', stderr);
  567. #endif /* def MPZ_TEST_DIV */
  568.  
  569.     (void)settupleitem(z, 0, (object *)x);
  570.     (void)settupleitem(z, 1, (object *)y);
  571.     
  572.     return z;
  573. } /* mpz_div_and_mod() */
  574.  
  575. static object *
  576. mpz_power(a, b, m)
  577.     mpzobject *a;
  578.     mpzobject *b;
  579.         mpzobject *m;
  580. {
  581.     mpzobject *z;
  582.     int cmpres;
  583.     long int longtmp1, longtmp2;
  584.  
  585.      if ((object *)m!=Py_None)
  586.        {
  587.          mpzobject *z2;
  588.         INCREF(Py_None);
  589.          z=(mpzobject *)mpz_power(a, b, (mpzobject *)Py_None);
  590.         DECREF(Py_None);
  591.          if (z==NULL) return((object *)z);
  592.          z2=(mpzobject *)mpz_remainder(z, m);
  593.          DECREF(z);
  594.          return((object *)z2);
  595.        }        
  596.  
  597.     if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
  598.         /* the gnu-mp lib sets pow(0,0) to 0, we to 1 */
  599.  
  600.         INCREF(mpz_value_one);
  601.         return (object *)mpz_value_one;
  602.     }
  603.         
  604.     if (cmpres < 0) {
  605.         err_setstr(ValueError, "mpz.pow to negative exponent");
  606.         return NULL;
  607.     }
  608.  
  609.     if ((cmpres = mpz_cmp_ui(&a->mpz, (unsigned long int)0)) == 0) {
  610.         /* the base is 0 */
  611.  
  612.         INCREF(mpz_value_zero);
  613.         return (object *)mpz_value_zero;
  614.     }
  615.     else if (cmpres > 0
  616.          && mpz_cmp_ui(&a->mpz, (unsigned long int)1) == 0) {
  617.         /* the base is 1 */
  618.  
  619.         INCREF(mpz_value_one);
  620.         return (object *)mpz_value_one;
  621.     }
  622.     else if (cmpres < 0
  623.          && mpz_cmp_si(&a->mpz, (long int)-1) == 0) {
  624.  
  625.         MP_INT tmpmpz;
  626.         /* the base is -1: pow(-1, any) == 1,-1 for even,uneven b */
  627.         /* XXX this code needs to be optimized: what's better?
  628.            mpz_mmod_ui or mpz_mod_2exp, I choose for the latter
  629.            for *un*obvious reasons */
  630.  
  631.         /* is the exponent even? */
  632.         mpz_init(&tmpmpz);
  633.  
  634.         /* look to the remainder after a division by (1 << 1) */
  635.         mpz_mod_2exp(&tmpmpz, &b->mpz, (unsigned long int)1);
  636.  
  637.         if (mpz_cmp_ui(&tmpmpz, (unsigned int)0) == 0) {
  638.             mpz_clear(&tmpmpz);
  639.             INCREF(mpz_value_one);
  640.             return (object *)mpz_value_one;
  641.         }
  642.         mpz_clear(&tmpmpz);
  643.         INCREF(mpz_value_mone);
  644.         return (object *)mpz_value_mone;
  645.     }
  646.  
  647. #ifdef MPZ_LIB_DOES_CHECKING
  648.     /* check if it's doable: sizeof(exp) > sizeof(long) &&
  649.        abs(base) > 1 ?? --> No Way */
  650.     if (mpz_size(&b->mpz) > 1)
  651.         return (object *)err_nomem();
  652. #else /* def MPZ_LIB_DOES_CHECKING */
  653.     /* wet finger method */
  654.     if (mpz_cmp_ui(&b->mpz, (unsigned long int)0x10000) >= 0) {
  655.         err_setstr(ValueError, "mpz.pow outrageous exponent");
  656.         return NULL;
  657.     }
  658. #endif /* def MPZ_LIB_DOES_CHECKING else */
  659.  
  660.     if ((z = newmpzobject()) == NULL)
  661.         return NULL;
  662.     
  663.     mpz_pow_ui(&z->mpz, &a->mpz, mpz_get_ui(&b->mpz));
  664.     
  665.     return (object *)z;
  666. } /* mpz_power() */
  667.  
  668.  
  669. static object *
  670. mpz_negative(v)
  671.     mpzobject *v;
  672. {
  673.     mpzobject *z;
  674.  
  675.     
  676. #ifdef MPZ_SPARE_MALLOC
  677.     if (mpz_cmp_ui(&v->mpz, (unsigned long int)0) == 0) {
  678.         /* -0 == 0 */
  679.         INCREF(v);
  680.         return (object *)v;
  681.     }
  682. #endif /* def MPZ_SPARE_MALLOC */
  683.  
  684.     if ((z = newmpzobject()) == NULL)
  685.         return NULL;
  686.  
  687.     mpz_neg(&z->mpz, &v->mpz);
  688.     return (object *)z;
  689. } /* mpz_negative() */
  690.  
  691.  
  692. static object *
  693. mpz_positive(v)
  694.     mpzobject *v;
  695. {
  696.     INCREF(v);
  697.     return (object *)v;
  698. } /* mpz_positive() */
  699.  
  700.  
  701. static object *
  702. mpz_absolute(v)
  703.     mpzobject *v;
  704. {
  705.     mpzobject *z;
  706.  
  707.     
  708.     if (mpz_cmp_ui(&v->mpz, (unsigned long int)0) >= 0) {
  709.         INCREF(v);
  710.         return (object *)v;
  711.     }
  712.  
  713.     if ((z = newmpzobject()) == NULL)
  714.         return NULL;
  715.  
  716.     mpz_neg(&z->mpz, &v->mpz);
  717.     return (object *)z;
  718. } /* mpz_absolute() */
  719.  
  720. static int
  721. mpz_nonzero(v)
  722.     mpzobject *v;
  723. {
  724.     return mpz_cmp_ui(&v->mpz, (unsigned long int)0) != 0;
  725. } /* mpz_nonzero() */
  726.         
  727. static object *
  728. mpz_invert(v)
  729.     mpzobject *v;
  730. {
  731.     mpzobject *z;
  732.  
  733.  
  734.     /* I think mpz_com does exactly what needed */
  735.     if ((z = newmpzobject()) == NULL)
  736.         return NULL;
  737.  
  738.     mpz_com(&z->mpz, &v->mpz);
  739.     return (object *)z;
  740. } /* mpz_invert() */
  741.  
  742. static object *
  743. mpz_lshift(a, b)
  744.     mpzobject *a;
  745.     mpzobject *b;
  746. {
  747.     int cmpres;
  748.     mpzobject *z;
  749.  
  750.  
  751.     if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
  752.         /* a << 0 == a */
  753.         INCREF(a);
  754.         return (object *)a;
  755.     }
  756.  
  757.     if (cmpres < 0) {
  758.         err_setstr(ValueError, "mpz.<< negative shift count");
  759.         return NULL;
  760.     }
  761.  
  762. #ifdef MPZ_LIB_DOES_CHECKING
  763.     if (mpz_size(&b->mpz) > 1)
  764.         return (object *)err_nomem();
  765. #else /* def MPZ_LIB_DOES_CHECKING */
  766.     /* wet finger method */
  767.     if (mpz_cmp_ui(&b->mpz, (unsigned long int)0x10000) >= 0) {
  768.         err_setstr(ValueError, "mpz.<< outrageous shift count");
  769.         return NULL;
  770.     }
  771. #endif /* def MPZ_LIB_DOES_CHECKING else */
  772.  
  773.     if ((z = newmpzobject()) == NULL)
  774.         return NULL;
  775.  
  776.     mpz_mul_2exp(&z->mpz, &a->mpz, mpz_get_ui(&b->mpz));
  777.     return (object *)z;
  778. } /* mpz_lshift() */
  779.  
  780. static object *
  781. mpz_rshift(a, b)
  782.     mpzobject *a;
  783.     mpzobject *b;
  784. {
  785.     int cmpres;
  786.     mpzobject *z;
  787.  
  788.  
  789.     if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
  790.         /* a >> 0 == a */
  791.         INCREF(a);
  792.         return (object *)a;
  793.     }
  794.  
  795.     if (cmpres < 0) {
  796.         err_setstr(ValueError, "mpz.>> negative shift count");
  797.         return NULL;
  798.     }
  799.  
  800.     if (mpz_size(&b->mpz) > 1)
  801.         return (object *)err_nomem();
  802.  
  803.     if ((z = newmpzobject()) == NULL)
  804.         return NULL;
  805.  
  806.     mpz_div_2exp(&z->mpz, &a->mpz, mpz_get_ui(&b->mpz));
  807.     return (object *)z;
  808. } /* mpz_rshift() */
  809.  
  810. static object *
  811. mpz_andfunc(a, b)
  812.     mpzobject *a;
  813.     mpzobject *b;
  814. {
  815.     mpzobject *z;
  816.  
  817.  
  818.     if ((z = newmpzobject()) == NULL)
  819.         return NULL;
  820.  
  821.     mpz_and(&z->mpz, &a->mpz, &b->mpz);
  822.     return (object *)z;
  823. } /* mpz_andfunc() */
  824.  
  825. /* hack Hack HAck HACk HACK, XXX this code is dead slow */
  826. void
  827. mpz_xor(res, op1, op2)
  828.     MP_INT *res;
  829.     const MP_INT *op1;
  830.     const MP_INT *op2;
  831. {
  832.     MP_INT tmpmpz;
  833.     
  834.     mpz_init(&tmpmpz);
  835.  
  836.     mpz_and(res, op1, op2);
  837.     mpz_com(&tmpmpz, res);
  838.     mpz_ior(res, op1, op2);
  839.     mpz_and(res, res, &tmpmpz);
  840.  
  841.     mpz_clear(&tmpmpz);
  842. } /* mpz_xor() HACK */
  843.  
  844. static object *
  845. mpz_xorfunc(a, b)
  846.     mpzobject *a;
  847.     mpzobject *b;
  848. {
  849.     mpzobject *z;
  850.  
  851.  
  852.     if ((z = newmpzobject()) == NULL)
  853.         return NULL;
  854.  
  855.     mpz_xor(&z->mpz, &a->mpz, &b->mpz);
  856.     return (object *)z;
  857. } /* mpz_xorfunc() */
  858.  
  859. static object *
  860. mpz_orfunc(a, b)
  861.     mpzobject *a;
  862.     mpzobject *b;
  863. {
  864.     mpzobject *z;
  865.  
  866.  
  867.     if ((z = newmpzobject()) == NULL)
  868.         return NULL;
  869.  
  870.     mpz_ior(&z->mpz, &a->mpz, &b->mpz);
  871.     return (object *)z;
  872. } /* mpz_orfunc() */
  873.  
  874. /* MPZ initialisation */
  875.  
  876. #include "longintrepr.h"
  877.  
  878. static object *
  879. MPZ_mpz(self, args)
  880.     object *self;
  881.     object *args;
  882. {
  883.     mpzobject *mpzp;
  884.     object *objp;
  885.  
  886.  
  887. #ifdef MPZ_DEBUG
  888.     fputs("MPZ_mpz() called...\n", stderr);
  889. #endif /* def MPZ_DEBUG */
  890.  
  891.     if (!getargs(args, "O", &objp))
  892.         return NULL;
  893.  
  894.     /* at least we know it's some object */
  895.     /* note DON't DECREF args NEITHER objp */
  896.  
  897.     if (is_intobject(objp)) {
  898.         long lval;
  899.  
  900.         if (!getargs(objp, "l", &lval))
  901.             return NULL;
  902.         
  903.         if (lval == (long)0) {
  904.             INCREF(mpz_value_zero);
  905.             mpzp = mpz_value_zero;
  906.         }
  907.         else if (lval == (long)1) {
  908.             INCREF(mpz_value_one);
  909.             mpzp = mpz_value_one;
  910.         }            
  911.         else if ((mpzp = newmpzobject()) == NULL)
  912.             return NULL;
  913.         else mpz_set_si(&mpzp->mpz, lval);
  914.     }
  915.     else if (is_longobject(objp)) {
  916.         MP_INT mplongdigit;
  917.         int i;
  918.         unsigned char isnegative;
  919.         
  920.  
  921.         if ((mpzp = newmpzobject()) == NULL)
  922.             return NULL;
  923.  
  924.         mpz_set_si(&mpzp->mpz, 0L);
  925.         mpz_init(&mplongdigit);
  926.         
  927.         /* how we're gonna handle this? */
  928.         if (isnegative = ((i = ((longobject *)objp)->ob_size) < 0) )
  929.             i = -i;
  930.  
  931.         while (i--) {
  932.             mpz_set_ui(&mplongdigit,
  933.                    (unsigned long)
  934.                    ((longobject *)objp)->ob_digit[i]);
  935.             mpz_mul_2exp(&mplongdigit,&mplongdigit,
  936.                      (unsigned long int)i * SHIFT);
  937.             mpz_ior(&mpzp->mpz, &mpzp->mpz, &mplongdigit);
  938.         }
  939.  
  940.         if (isnegative)
  941.             mpz_neg(&mpzp->mpz, &mpzp->mpz);
  942.  
  943.         /* get rid of allocation for tmp variable */
  944.         mpz_clear(&mplongdigit);
  945.     }
  946.     else if (is_stringobject(objp)) {
  947.         char *cp;
  948.         int len;
  949.         MP_INT mplongdigit;
  950.         
  951.         if (!getargs(objp, "s#", &cp, &len))
  952.             return NULL;
  953.  
  954.         if ((mpzp = newmpzobject()) == NULL)
  955.             return NULL;
  956.  
  957.         mpz_set_si(&mpzp->mpz, 0L);
  958.         mpz_init(&mplongdigit);
  959.         
  960.         /* let's do it the same way as with the long conversion:
  961.            without thinking how it can be faster (-: :-) */
  962.  
  963.         cp += len;
  964.         while (len--) {
  965.             mpz_set_ui(&mplongdigit, (unsigned long)*--cp );
  966.             mpz_mul_2exp(&mplongdigit,&mplongdigit,
  967.                      (unsigned long int)len * 8);
  968.             mpz_ior(&mpzp->mpz, &mpzp->mpz, &mplongdigit);
  969.         }
  970.  
  971.         /* get rid of allocation for tmp variable */
  972.         mpz_clear(&mplongdigit);
  973.     }
  974.     else if (is_mpzobject(objp)) {
  975.         INCREF(objp);
  976.         mpzp = (mpzobject *)objp;
  977.     }
  978.     else {
  979.         err_setstr(TypeError,
  980.                "mpz.mpz() expects integer, long, string or mpz object argument");
  981.         return NULL;
  982.     }
  983.  
  984.  
  985. #ifdef MPZ_DEBUG
  986.     fputs("MPZ_mpz: created mpz=", stderr);
  987.     mpz_out_str(stderr, 10, &mpzp->mpz);
  988.     putc('\n', stderr);
  989. #endif /* def MPZ_DEBUG */
  990.     return (object *)mpzp;
  991. } /* MPZ_mpz() */
  992.  
  993. static mpzobject *
  994. mpz_mpzcoerce(z)
  995.     object *z;
  996. {
  997.     /* shortcut: 9 out of 10 times the type is already ok */
  998.     if (is_mpzobject(z)) {
  999.         INCREF(z);
  1000.         return (mpzobject *)z;    /* coercion succeeded */
  1001.     }
  1002.  
  1003.     /* what types do we accept?: intobjects and longobjects */
  1004.     if (is_intobject(z) || is_longobject(z))
  1005.         return (mpzobject *)MPZ_mpz((object *)NULL, z);
  1006.  
  1007.     err_setstr(TypeError, "number coercion (to mpzobject) failed");
  1008.     return NULL;
  1009. } /* mpz_mpzcoerce() */
  1010.     
  1011. /* Forward */
  1012. static void mpz_divm PROTO((MP_INT *res, const MP_INT *num,
  1013.                 const MP_INT *den, const MP_INT *mod));
  1014.  
  1015. static object *
  1016. MPZ_powm(self, args)
  1017.     object *self;
  1018.     object *args;
  1019. {
  1020.     object *base, *exp, *mod;
  1021.     mpzobject *mpzbase = NULL, *mpzexp = NULL, *mpzmod = NULL;
  1022.     mpzobject *z;
  1023.     int tstres;
  1024.  
  1025.     
  1026.     if (!getargs(args, "(OOO)", &base, &exp, &mod))
  1027.         return NULL;
  1028.  
  1029.     if ((mpzbase = mpz_mpzcoerce(base)) == NULL
  1030.         || (mpzexp = mpz_mpzcoerce(exp)) == NULL
  1031.         || (mpzmod = mpz_mpzcoerce(mod)) == NULL
  1032.         || (z = newmpzobject()) == NULL) {
  1033.         XDECREF(mpzbase);
  1034.         XDECREF(mpzexp);
  1035.         XDECREF(mpzmod);
  1036.         return NULL;
  1037.     }
  1038.  
  1039.     if ((tstres=mpz_cmp_ui(&mpzexp->mpz, (unsigned long int)0)) == 0) {
  1040.         INCREF(mpz_value_one);
  1041.         return (object *)mpz_value_one;
  1042.     }
  1043.  
  1044.     if (tstres < 0) {
  1045.         MP_INT absexp;
  1046.         /* negative exp */
  1047.  
  1048.         mpz_init_set(&absexp, &mpzexp->mpz);
  1049.         mpz_abs(&absexp, &absexp);
  1050.         mpz_powm(&z->mpz, &mpzbase->mpz, &absexp, &mpzmod->mpz);
  1051.  
  1052.         mpz_divm(&z->mpz, &mpz_value_one->mpz, &z->mpz, &mpzmod->mpz);
  1053.         
  1054.         mpz_clear(&absexp);
  1055.     }
  1056.     else {
  1057.         mpz_powm(&z->mpz, &mpzbase->mpz, &mpzexp->mpz, &mpzmod->mpz);
  1058.     }
  1059.         
  1060.     DECREF(mpzbase);
  1061.     DECREF(mpzexp);
  1062.     DECREF(mpzmod);
  1063.  
  1064.     return (object *)z;
  1065. } /* MPZ_powm() */
  1066.  
  1067.  
  1068. static object *
  1069. MPZ_gcd(self, args)
  1070.     object *self;
  1071.     object *args;
  1072. {
  1073.     object *op1, *op2;
  1074.     mpzobject *mpzop1 = NULL, *mpzop2 = NULL;
  1075.     mpzobject *z;
  1076.  
  1077.     
  1078.     if (!getargs(args, "(OO)", &op1, &op2))
  1079.         return NULL;
  1080.  
  1081.     if ((mpzop1 = mpz_mpzcoerce(op1)) == NULL
  1082.         || (mpzop2 = mpz_mpzcoerce(op2)) == NULL
  1083.         || (z = newmpzobject()) == NULL) {
  1084.         XDECREF(mpzop1);
  1085.         XDECREF(mpzop2);
  1086.         return NULL;
  1087.     }
  1088.  
  1089.     /* ok, we have three mpzobjects, and an initialised result holder */
  1090.     mpz_gcd(&z->mpz, &mpzop1->mpz, &mpzop2->mpz);
  1091.  
  1092.     DECREF(mpzop1);
  1093.     DECREF(mpzop2);
  1094.  
  1095.     return (object *)z;
  1096. } /* MPZ_gcd() */
  1097.  
  1098.  
  1099. static object *
  1100. MPZ_gcdext(self, args)
  1101.     object *self;
  1102.     object *args;
  1103. {
  1104.     object *op1, *op2, *z = NULL;
  1105.     mpzobject *mpzop1 = NULL, *mpzop2 = NULL;
  1106.     mpzobject *g = NULL, *s = NULL, *t = NULL;
  1107.  
  1108.     
  1109.     if (!getargs(args, "(OO)", &op1, &op2))
  1110.         return NULL;
  1111.  
  1112.     if ((mpzop1 = mpz_mpzcoerce(op1)) == NULL
  1113.         || (mpzop2 = mpz_mpzcoerce(op2)) == NULL
  1114.         || (z = newtupleobject(3)) == NULL
  1115.         || (g = newmpzobject()) == NULL
  1116.         || (s = newmpzobject()) == NULL
  1117.         || (t = newmpzobject()) == NULL) {
  1118.         XDECREF(mpzop1);
  1119.         XDECREF(mpzop2);
  1120.         XDECREF(z);
  1121.         XDECREF(g);
  1122.         XDECREF(s);
  1123.         /*XDECREF(t);*/
  1124.         return NULL;
  1125.     }
  1126.  
  1127.     mpz_gcdext(&g->mpz, &s->mpz, &t->mpz, &mpzop1->mpz, &mpzop2->mpz);
  1128.  
  1129.     DECREF(mpzop1);
  1130.     DECREF(mpzop2);
  1131.  
  1132.     (void)settupleitem(z, 0, (object *)g);
  1133.     (void)settupleitem(z, 1, (object *)s);
  1134.     (void)settupleitem(z, 2, (object *)t);
  1135.  
  1136.     return (object *)z;
  1137. } /* MPZ_gcdext() */
  1138.  
  1139.  
  1140. static object *
  1141. MPZ_sqrt(self, args)
  1142.     object *self;
  1143.     object *args;
  1144. {
  1145.     object *op;
  1146.     mpzobject *mpzop = NULL;
  1147.     mpzobject *z;
  1148.  
  1149.     
  1150.     if (!getargs(args, "O", &op))
  1151.         return NULL;
  1152.  
  1153.     if ((mpzop = mpz_mpzcoerce(op)) == NULL
  1154.         || (z = newmpzobject()) == NULL) {
  1155.         XDECREF(mpzop);
  1156.         return NULL;
  1157.     }
  1158.  
  1159.     mpz_sqrt(&z->mpz, &mpzop->mpz);
  1160.  
  1161.     DECREF(mpzop);
  1162.  
  1163.     return (object *)z;
  1164. } /* MPZ_sqrt() */
  1165.  
  1166.  
  1167. static object *
  1168. MPZ_sqrtrem(self, args)
  1169.     object *self;
  1170.     object *args;
  1171. {
  1172.     object *op, *z = NULL;
  1173.     mpzobject *mpzop = NULL;
  1174.     mpzobject *root = NULL, *rem = NULL;
  1175.  
  1176.     
  1177.     if (!getargs(args, "O", &op))
  1178.         return NULL;
  1179.  
  1180.     if ((mpzop = mpz_mpzcoerce(op)) == NULL
  1181.         || (z = newtupleobject(2)) == NULL
  1182.         || (root = newmpzobject()) == NULL
  1183.         || (rem = newmpzobject()) == NULL) {
  1184.         XDECREF(mpzop);
  1185.         XDECREF(z);
  1186.         XDECREF(root);
  1187.         /*XDECREF(rem);*/
  1188.         return NULL;
  1189.     }
  1190.  
  1191.     mpz_sqrtrem(&root->mpz, &rem->mpz, &mpzop->mpz);
  1192.  
  1193.     DECREF(mpzop);
  1194.  
  1195.     (void)settupleitem(z, 0, (object *)root);
  1196.     (void)settupleitem(z, 1, (object *)rem);
  1197.  
  1198.     return (object *)z;
  1199. } /* MPZ_sqrtrem() */
  1200.  
  1201.  
  1202. static void
  1203. #if __STDC__
  1204. mpz_divm(MP_INT *res, const MP_INT *num, const MP_INT *den, const MP_INT *mod)
  1205. #else
  1206. mpz_divm(res, num, den, mod)
  1207.     MP_INT *res;
  1208.     const MP_INT *num;
  1209.     const MP_INT *den;
  1210.     const MP_INT *mod;
  1211. #endif
  1212. {
  1213.     MP_INT s0, s1, q, r, x, d0, d1;
  1214.  
  1215.     mpz_init_set(&s0, num);
  1216.     mpz_init_set_ui(&s1, 0);
  1217.     mpz_init(&q);
  1218.     mpz_init(&r);
  1219.     mpz_init(&x);
  1220.     mpz_init_set(&d0, den);
  1221.     mpz_init_set(&d1, mod);
  1222.  
  1223.     while (d1.size != 0) {
  1224.         mpz_divmod(&q, &r, &d0, &d1);
  1225.         mpz_set(&d0, &d1);
  1226.         mpz_set(&d1, &r);
  1227.  
  1228.         mpz_mul(&x, &s1, &q);
  1229.         mpz_sub(&x, &s0, &x);
  1230.         mpz_set(&s0, &s1);
  1231.         mpz_set(&s1, &x);
  1232.     }
  1233.  
  1234.     if (d0.size != 1 || d0.d[0] != 1)
  1235.         res->size = 0;    /* trouble: the gcd != 1; set s to zero */
  1236.     else {
  1237. #ifdef MPZ_MDIV_BUG
  1238.         /* watch out here! first check the signs, and then perform
  1239.            the mpz_mod() since mod could point to res */
  1240.         if ((s0.size < 0) != (mod->size < 0)) {
  1241.             mpz_mod(res, &s0, mod);
  1242.  
  1243.             if (res->size)
  1244.                 mpz_add(res, res, mod);
  1245.         }
  1246.         else
  1247.             mpz_mod(res, &s0, mod);
  1248.         
  1249. #else /* def MPZ_MDIV_BUG */
  1250.         mpz_mmod(res, &s0, mod);
  1251. #endif /* def MPZ_MDIV_BUG else */
  1252.     }
  1253.  
  1254.     mpz_clear(&s0);
  1255.     mpz_clear(&s1);
  1256.     mpz_clear(&q);
  1257.     mpz_clear(&r);
  1258.     mpz_clear(&x);
  1259.     mpz_clear(&d0);
  1260.     mpz_clear(&d1);
  1261. } /* mpz_divm() */
  1262.  
  1263.  
  1264. static object *
  1265. MPZ_divm(self, args)
  1266.     object *self;
  1267.     object *args;
  1268. {
  1269.     object *num, *den, *mod;
  1270.     mpzobject *mpznum, *mpzden, *mpzmod = NULL;
  1271.     mpzobject *z = NULL;
  1272.  
  1273.     
  1274.     if (!getargs(args, "(OOO)", &num, &den, &mod))
  1275.         return NULL;
  1276.  
  1277.     if ((mpznum = mpz_mpzcoerce(num)) == NULL
  1278.         || (mpzden = mpz_mpzcoerce(den)) == NULL
  1279.         || (mpzmod = mpz_mpzcoerce(mod)) == NULL
  1280.         || (z = newmpzobject()) == NULL ) {
  1281.         XDECREF(mpznum);
  1282.         XDECREF(mpzden);
  1283.         XDECREF(mpzmod);
  1284.         return NULL;
  1285.     }
  1286.     
  1287.     mpz_divm(&z->mpz, &mpznum->mpz, &mpzden->mpz, &mpzmod->mpz);
  1288.  
  1289.     DECREF(mpznum);
  1290.     DECREF(mpzden);
  1291.     DECREF(mpzmod);
  1292.  
  1293.     if (mpz_cmp_ui(&z->mpz, (unsigned long int)0) == 0) {
  1294.         DECREF(z);
  1295.         err_setstr(ValueError, "gcd(den, mod) != 1 or num == 0");
  1296.         return NULL;
  1297.     }
  1298.  
  1299.     return (object *)z;
  1300. } /* MPZ_divm() */
  1301.  
  1302.  
  1303. /* MPZ methods-as-attributes */
  1304. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1305. static object *
  1306. mpz_int(self, args)
  1307.     mpzobject *self;
  1308.     object *args;
  1309. #else /* def MPZ_CONVERSIONS_AS_METHODS */
  1310. static object *
  1311. mpz_int(self)
  1312.     mpzobject *self;
  1313. #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
  1314. {
  1315.     long sli;
  1316.  
  1317.  
  1318. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1319.     if (!getnoarg(args))
  1320.         return NULL;
  1321. #endif /* def MPZ_CONVERSIONS_AS_METHODS */
  1322.  
  1323.     if (mpz_size(&self->mpz) > 1
  1324.         || (sli = (long)mpz_get_ui(&self->mpz)) < (long)0 ) {
  1325.         err_setstr(ValueError, "mpz.int() arg too long to convert");
  1326.         return NULL;
  1327.     }
  1328.  
  1329.     if (mpz_cmp_ui(&self->mpz, (unsigned long)0) < 0)
  1330.         sli = -sli;
  1331.  
  1332.     return newintobject(sli);
  1333. } /* mpz_int() */
  1334.     
  1335. static object *
  1336. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1337. mpz_long(self, args)
  1338.     mpzobject *self;
  1339.     object *args;
  1340. #else /* def MPZ_CONVERSIONS_AS_METHODS */
  1341. mpz_long(self)
  1342.     mpzobject *self;
  1343. #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
  1344. {
  1345.     int i, isnegative;
  1346.     unsigned long int uli;
  1347.     longobject *longobjp;
  1348.     int ldcount;
  1349.     int bitpointer, newbitpointer;
  1350.     MP_INT mpzscratch;
  1351.  
  1352.  
  1353. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1354.     if (!getnoarg(args))
  1355.         return NULL;
  1356. #endif /* def MPZ_CONVERSIONS_AS_METHODS */
  1357.  
  1358.     /* determine length of python-long to be allocated */
  1359.     if ((longobjp = alloclongobject(i = (int)
  1360.                 ((mpz_size(&self->mpz) * BITS_PER_MP_LIMB
  1361.                   + SHIFT - 1) /
  1362.                  SHIFT))) == NULL)
  1363.         return NULL;
  1364.  
  1365.     /* determine sign, and copy self to scratch var */
  1366.     mpz_init_set(&mpzscratch, &self->mpz);
  1367.     if (isnegative = (mpz_cmp_ui(&self->mpz, (unsigned long int)0) < 0))
  1368.         mpz_neg(&mpzscratch, &mpzscratch);
  1369.  
  1370.     /* let those bits come, let those bits go,
  1371.        e.g. dismantle mpzscratch, build longobject */
  1372.  
  1373.     bitpointer = 0;        /* the number of valid bits in stock */
  1374.     newbitpointer = 0;
  1375.     ldcount = 0;        /* the python-long limb counter */
  1376.     uli = (unsigned long int)0;
  1377.     while (i--) {
  1378.         longobjp->ob_digit[ldcount] = uli & MASK;
  1379.  
  1380.         /* check if we've had enough bits for this digit */
  1381.         if (bitpointer < SHIFT) {
  1382.             uli = mpz_get_ui(&mpzscratch);
  1383.             longobjp->ob_digit[ldcount] |=
  1384.                 (uli << bitpointer) & MASK;
  1385.             uli >>= SHIFT-bitpointer;
  1386.             bitpointer += BITS_PER_MP_LIMB;
  1387.             mpz_div_2exp(&mpzscratch, &mpzscratch,
  1388.                      BITS_PER_MP_LIMB);
  1389.         }
  1390.         else
  1391.             uli >>= SHIFT;
  1392.         bitpointer -= SHIFT;
  1393.         ldcount++;
  1394.     }
  1395.  
  1396.     assert(mpz_cmp_ui(&mpzscratch, (unsigned long int)0) == 0);
  1397.     mpz_clear(&mpzscratch);
  1398.     assert(ldcount <= longobjp->ob_size);
  1399.  
  1400.     /* long_normalize() is file-static */
  1401.     /* longobjp = long_normalize(longobjp); */
  1402.     while (ldcount > 0 && longobjp->ob_digit[ldcount-1] == 0)
  1403.         ldcount--;
  1404.     longobjp->ob_size = ldcount;
  1405.     
  1406.  
  1407.     if (isnegative)
  1408.         longobjp->ob_size = -longobjp->ob_size;
  1409.  
  1410.     return (object *)longobjp;
  1411.     
  1412. } /* mpz_long() */
  1413.  
  1414.  
  1415. /* I would have avoided pow() anyways, so ... */
  1416. static const double multiplier = 256.0 * 256.0 * 256.0 * 256.0;
  1417.     
  1418. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1419. static object *
  1420. mpz_float(self, args)
  1421.     mpzobject *self;
  1422.     object *args;
  1423. #else /* def MPZ_CONVERSIONS_AS_METHODS */
  1424. static object *
  1425. mpz_float(self)
  1426.     mpzobject *self;
  1427. #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
  1428. {
  1429.     int i, isnegative;
  1430.     double x;
  1431.     double mulstate;
  1432.     MP_INT mpzscratch;
  1433.  
  1434.  
  1435. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1436.     if (!getnoarg(args))
  1437.         return NULL;
  1438. #endif /* def MPZ_CONVERSIONS_AS_METHODS */
  1439.  
  1440.     i = (int)mpz_size(&self->mpz);
  1441.     
  1442.     /* determine sign, and copy abs(self) to scratch var */
  1443.     if (isnegative = (mpz_cmp_ui(&self->mpz, (unsigned long int)0) < 0)) {
  1444.         mpz_init(&mpzscratch);
  1445.         mpz_neg(&mpzscratch, &self->mpz);
  1446.     }
  1447.     else
  1448.         mpz_init_set(&mpzscratch, &self->mpz);
  1449.  
  1450.     /* let those bits come, let those bits go,
  1451.        e.g. dismantle mpzscratch, build floatobject */
  1452.  
  1453.     x = 0.0;
  1454.     mulstate = 1.0;
  1455.     while (i--) {
  1456.         x += mulstate * mpz_get_ui(&mpzscratch);
  1457.         mulstate *= multiplier;
  1458.         mpz_div_2exp(&mpzscratch, &mpzscratch, BITS_PER_MP_LIMB);
  1459.     }
  1460.  
  1461.     assert(mpz_cmp_ui(&mpzscratch, (unsigned long int)0) == 0);
  1462.     mpz_clear(&mpzscratch);
  1463.  
  1464.     if (isnegative)
  1465.         x = -x;
  1466.  
  1467.     return newfloatobject(x);
  1468.     
  1469. } /* mpz_float() */
  1470.  
  1471. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1472. static object *
  1473. mpz_hex(self, args)
  1474.     mpzobject *self;
  1475.     object *args;
  1476. #else /* def MPZ_CONVERSIONS_AS_METHODS */
  1477. static object *
  1478. mpz_hex(self)
  1479.     mpzobject *self;
  1480. #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
  1481. {
  1482. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1483.     if (!getnoarg(args))
  1484.         return NULL;
  1485. #endif /* def MPZ_CONVERSIONS_AS_METHODS */
  1486.     
  1487.     return mpz_format(self, 16, (unsigned char)1);
  1488. } /* mpz_hex() */
  1489.     
  1490. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1491. static object *
  1492. mpz_oct(self, args)
  1493.     mpzobject *self;
  1494.     object *args;
  1495. #else /* def MPZ_CONVERSIONS_AS_METHODS */
  1496. static object *
  1497. mpz_oct(self)
  1498.     mpzobject *self;
  1499. #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
  1500. {
  1501. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1502.     if (!getnoarg(args))
  1503.         return NULL;
  1504. #endif /* def MPZ_CONVERSIONS_AS_METHODS */
  1505.     
  1506.     return mpz_format(self, 8, (unsigned char)1);
  1507. } /* mpz_oct() */
  1508.     
  1509. static object *
  1510. mpz_binary(self, args)
  1511.     mpzobject *self;
  1512.     object *args;
  1513. {
  1514.     int size;
  1515.     stringobject *strobjp;
  1516.     char *cp;
  1517.     MP_INT mp;
  1518.     unsigned long ldigit;
  1519.     
  1520.     if (!getnoarg(args))
  1521.         return NULL;
  1522.  
  1523.     if (mpz_cmp_ui(&self->mpz, (unsigned long int)0) < 0) {
  1524.         err_setstr(ValueError, "mpz.binary() arg must be >= 0");
  1525.         return NULL;
  1526.     }
  1527.  
  1528.     mpz_init_set(&mp, &self->mpz);
  1529.     size = (int)mpz_size(&mp);
  1530.  
  1531.     if ((strobjp = (stringobject *)
  1532.          newsizedstringobject((char *)0,
  1533.                   size * sizeof (unsigned long int))) == NULL)
  1534.         return NULL;
  1535.  
  1536.     /* get the beginning of the string memory and start copying things */
  1537.     cp = GETSTRINGVALUE(strobjp);
  1538.  
  1539.     /* this has been programmed using a (fairly) decent lib-i/f it could
  1540.        be must faster if we looked into the GMP lib */
  1541.     while (size--) {
  1542.         ldigit = mpz_get_ui(&mp);
  1543.         mpz_div_2exp(&mp, &mp, BITS_PER_MP_LIMB);
  1544.         *cp++ = (unsigned char)(ldigit & 0xFF);
  1545.         *cp++ = (unsigned char)((ldigit >>= 8) & 0xFF);
  1546.         *cp++ = (unsigned char)((ldigit >>= 8) & 0xFF);
  1547.         *cp++ = (unsigned char)((ldigit >>= 8) & 0xFF);
  1548.     }
  1549.  
  1550.     while (strobjp->ob_size && !*--cp)
  1551.         strobjp->ob_size--;
  1552.  
  1553.     return (object *)strobjp;
  1554. } /* mpz_binary() */
  1555.     
  1556.  
  1557. static struct methodlist mpz_methods[] = {
  1558. #ifdef MPZ_CONVERSIONS_AS_METHODS
  1559.     {"int",            mpz_int},
  1560.     {"long",        mpz_long},
  1561.     {"float",        mpz_float},
  1562.     {"hex",            mpz_hex},
  1563.     {"oct",            mpz_oct},
  1564. #endif /* def MPZ_CONVERSIONS_AS_METHODS */
  1565.     {"binary",        (object *(*)(object *, object *))mpz_binary},
  1566.     {NULL,            NULL}        /* sentinel */
  1567. };
  1568.  
  1569. static object *
  1570. mpz_getattr(self, name)
  1571.     mpzobject *self;
  1572.     char *name;
  1573. {
  1574.     return findmethod(mpz_methods, (object *)self, name);
  1575. } /* mpz_getattr() */
  1576.  
  1577.  
  1578. static int
  1579. mpz_coerce(pv, pw)
  1580.     object **pv;
  1581.     object **pw;
  1582. {
  1583.     object *z;
  1584.  
  1585. #ifdef MPZ_DEBUG
  1586.     fputs("mpz_coerce() called...\n", stderr);
  1587. #endif /* def MPZ_DEBUG */
  1588.  
  1589.     assert(is_mpzobject(*pv));
  1590.  
  1591.     /* always convert other arg to mpz value, except for floats */
  1592.     if (!is_floatobject(*pw)) {
  1593.         if ((z = (object *)mpz_mpzcoerce(*pw)) == NULL)
  1594.             return -1;    /* -1: an error always has been set */
  1595.         
  1596.         INCREF(*pv);
  1597.         *pw = z;
  1598.     }
  1599.     else {
  1600.         if ((z = mpz_float(*pv, NULL)) == NULL)
  1601.             return -1;
  1602.  
  1603.         INCREF(*pw);
  1604.         *pv = z;
  1605.     }
  1606.     return 0;        /* coercion succeeded */
  1607.  
  1608. } /* mpz_coerce() */
  1609.  
  1610.  
  1611. static object *
  1612. mpz_repr(v)
  1613.     object *v;
  1614. {
  1615.     return mpz_format(v, 10, (unsigned char)1);
  1616. } /* mpz_repr() */
  1617.  
  1618.  
  1619.  
  1620. #define UF (unaryfunc)
  1621. #define BF (binaryfunc)
  1622. #define TF (ternaryfunc)
  1623. #define IF (inquiry)
  1624. #define CF (coercion)
  1625.  
  1626. static number_methods mpz_as_number = {
  1627.     BF mpz_addition,    /*nb_add*/
  1628.     BF mpz_substract,    /*nb_subtract*/
  1629.     BF mpz_multiply,    /*nb_multiply*/
  1630.     BF mpz_divide,        /*nb_divide*/
  1631.     BF mpz_remainder,    /*nb_remainder*/
  1632.     BF mpz_div_and_mod,    /*nb_divmod*/
  1633.     TF mpz_power,        /*nb_power*/
  1634.     UF mpz_negative,    /*nb_negative*/
  1635.     UF mpz_positive,    /*tp_positive*/
  1636.     UF mpz_absolute,    /*tp_absolute*/
  1637.     IF mpz_nonzero,        /*tp_nonzero*/
  1638.     UF mpz_invert,        /*nb_invert*/
  1639.     BF mpz_lshift,        /*nb_lshift*/
  1640.     BF mpz_rshift,        /*nb_rshift*/
  1641.     BF mpz_andfunc,        /*nb_and*/
  1642.     BF mpz_xorfunc,        /*nb_xor*/
  1643.     BF mpz_orfunc,        /*nb_or*/
  1644.     CF mpz_coerce,        /*nb_coerce*/
  1645. #ifndef MPZ_CONVERSIONS_AS_METHODS
  1646.     UF mpz_int,        /*nb_int*/
  1647.     UF mpz_long,        /*nb_long*/
  1648.     UF mpz_float,        /*nb_float*/
  1649.     UF mpz_oct,        /*nb_oct*/
  1650.     UF mpz_hex,        /*nb_hex*/
  1651. #endif /* ndef MPZ_CONVERSIONS_AS_METHODS */
  1652. };
  1653.  
  1654. static typeobject MPZtype = {
  1655.     OB_HEAD_INIT(&Typetype)
  1656.     0,            /*ob_size*/
  1657.     "mpz",            /*tp_name*/
  1658.     sizeof(mpzobject),    /*tp_size*/
  1659.     0,            /*tp_itemsize*/
  1660.     /* methods */
  1661.     (destructor)mpz_dealloc, /*tp_dealloc*/
  1662.     0,            /*tp_print*/
  1663.     (getattrfunc)mpz_getattr, /*tp_getattr*/
  1664.     0,            /*tp_setattr*/
  1665.     (cmpfunc)mpz_compare,    /*tp_compare*/
  1666.     (reprfunc)mpz_repr,    /*tp_repr*/
  1667.         &mpz_as_number,     /*tp_as_number*/
  1668. };
  1669.  
  1670. /* List of functions exported by this module */
  1671.  
  1672. static struct methodlist mpz_functions[] = {
  1673. #if 0
  1674.     {initialiser_name,    MPZ_mpz},
  1675. #else /* 0 */
  1676.     /* until guido ``fixes'' struct methodlist */
  1677.     {(char *)initialiser_name,    MPZ_mpz},
  1678. #endif /* 0 else */    
  1679.     {"powm",        MPZ_powm},
  1680.     {"gcd",            MPZ_gcd},
  1681.     {"gcdext",        MPZ_gcdext},
  1682.     {"sqrt",        MPZ_sqrt},
  1683.     {"sqrtrem",        MPZ_sqrtrem},
  1684.     {"divm",        MPZ_divm},
  1685.     {NULL,            NULL}         /* Sentinel */
  1686. };
  1687.  
  1688.  
  1689. /* #define MP_TEST_ALLOC */
  1690.  
  1691. #ifdef MP_TEST_ALLOC
  1692. #define MP_TEST_SIZE        4
  1693. static const char mp_test_magic[MP_TEST_SIZE] = {'\xAA','\xAA','\xAA','\xAA'};
  1694. static mp_test_error( location )
  1695.     int *location;
  1696. {
  1697.     /* assumptions: *alloc returns address dividable by 4,
  1698.     mpz_* routines allocate in chunks dividable by four */
  1699.     fprintf(stderr, "MP_TEST_ERROR: location holds 0x%08d\n", *location );
  1700.     fatal("MP_TEST_ERROR");
  1701. } /* static mp_test_error() */
  1702. #define MP_EXTRA_ALLOC(size)    ((size) + MP_TEST_SIZE)
  1703. #define MP_SET_TEST(basep,size)    (void)memcpy( ((char *)(basep))+(size), mp_test_magic, MP_TEST_SIZE)
  1704. #define MP_DO_TEST(basep,size)    if ( !memcmp( ((char *)(basep))+(size), mp_test_magic, MP_TEST_SIZE ) ) \
  1705.                     ; \
  1706.                 else \
  1707.                     mp_test_error((int *)((char *)(basep) + size))
  1708. #else /* def MP_TEST_ALLOC */
  1709. #define MP_EXTRA_ALLOC(size)    (size)
  1710. #define MP_SET_TEST(basep,size)
  1711. #define MP_DO_TEST(basep,size)
  1712. #endif /* def MP_TEST_ALLOC else */
  1713.  
  1714. void *mp_allocate( alloc_size )
  1715.     size_t    alloc_size;
  1716. {
  1717.     void *res;
  1718.  
  1719. #ifdef MPZ_DEBUG
  1720.     fprintf(stderr, "mp_allocate  :                             size %ld\n",
  1721.         alloc_size);
  1722. #endif /* def MPZ_DEBUG */    
  1723.  
  1724.     if ( (res = malloc(MP_EXTRA_ALLOC(alloc_size))) == NULL )
  1725.         fatal("mp_allocate failure");
  1726.  
  1727. #ifdef MPZ_DEBUG
  1728.     fprintf(stderr, "mp_allocate  :     address 0x%08x\n", res);
  1729. #endif /* def MPZ_DEBUG */    
  1730.  
  1731.     MP_SET_TEST(res,alloc_size);
  1732.     
  1733.     return res;
  1734. } /* mp_allocate() */
  1735.  
  1736.  
  1737. void *mp_reallocate( ptr, old_size, new_size )
  1738.     void *ptr;
  1739.     size_t old_size;
  1740.     size_t new_size;
  1741. {
  1742.     void *res;
  1743.  
  1744. #ifdef MPZ_DEBUG
  1745.     fprintf(stderr, "mp_reallocate: old address 0x%08x, old size %ld\n",
  1746.         ptr, old_size);
  1747. #endif /* def MPZ_DEBUG */    
  1748.  
  1749.     MP_DO_TEST(ptr, old_size);
  1750.     
  1751.     if ( (res = realloc(ptr, MP_EXTRA_ALLOC(new_size))) == NULL )
  1752.         fatal("mp_reallocate failure");
  1753.  
  1754. #ifdef MPZ_DEBUG
  1755.     fprintf(stderr, "mp_reallocate: new address 0x%08x, new size %ld\n",
  1756.         res, new_size);
  1757. #endif /* def MPZ_DEBUG */    
  1758.  
  1759.     MP_SET_TEST(res, new_size);
  1760.  
  1761.     return res;
  1762. } /* mp_reallocate() */
  1763.  
  1764.  
  1765. void mp_free( ptr, size )
  1766.     void *ptr;
  1767.     size_t size;
  1768. {
  1769.  
  1770. #ifdef MPZ_DEBUG
  1771.     fprintf(stderr, "mp_free      : old address 0x%08x, old size %ld\n",
  1772.         ptr, size);
  1773. #endif /* def MPZ_DEBUG */    
  1774.  
  1775.     MP_DO_TEST(ptr, size);
  1776.     free(ptr);
  1777. } /* mp_free() */
  1778.  
  1779.  
  1780.  
  1781. /* Initialize this module. */
  1782.  
  1783. void
  1784. initmpz()
  1785. {
  1786. #ifdef MPZ_DEBUG
  1787.     fputs( "initmpz() called...\n", stderr );
  1788. #endif /* def MPZ_DEBUG */
  1789.  
  1790.     mp_set_memory_functions( mp_allocate, mp_reallocate, mp_free );
  1791.     (void)initmodule("mpz", mpz_functions);
  1792.  
  1793.     /* create some frequently used constants */
  1794.     if ((mpz_value_zero = newmpzobject()) == NULL)
  1795.         fatal("initmpz: can't initialize mpz constants");
  1796.     mpz_set_ui(&mpz_value_zero->mpz, (unsigned long int)0);
  1797.  
  1798.     if ((mpz_value_one = newmpzobject()) == NULL)
  1799.         fatal("initmpz: can't initialize mpz constants");
  1800.     mpz_set_ui(&mpz_value_one->mpz, (unsigned long int)1);
  1801.  
  1802.     if ((mpz_value_mone = newmpzobject()) == NULL)
  1803.         fatal("initmpz: can't initialize mpz constants");
  1804.     mpz_set_si(&mpz_value_mone->mpz, (long)-1);
  1805.  
  1806. } /* initmpz() */
  1807. #ifdef MAKEDUMMYINT
  1808. int _mpz_dummy_int;    /* XXX otherwise, we're .bss-less (DYNLOAD->Jack?) */
  1809. #endif /* def MAKEDUMMYINT */
  1810.